1 Task info and setup

This R Markdown script assesses models for the PAL (probabilistic associative learning) task of the EMBA project.

1.1 Some general settings

# number of simulations
nsim = 250

# set the seed
set.seed(2468)

1.2 Package versions

The following packages are used in this RMarkdown file:

## [1] "R version 4.5.1 (2025-06-13)"
## [1] "knitr version 1.50"
## [1] "ggplot2 version 3.5.2"
## [1] "brms version 2.22.0"
## [1] "bridgesampling version 1.1.2"
## [1] "tidyverse version 2.0.0"
## [1] "ggpubr version 0.6.1"
## [1] "ggrain version 0.0.4"
## [1] "bayesplot version 1.13.0"
## [1] "SBC version 0.3.0.9000"
## [1] "rstatix version 0.7.2"
## [1] "effectsize version 1.0.1"
## [1] "BayesFactor version 0.9.12.4.7"
## [1] "bayestestR version 0.17.0"

1.3 Evaluation guidelines

We perform prior predictive checks as proposed in Schad, Betancourt and Vasishth (2020) using the SBC package. To do so, we create simulated datasets where parameters are simulated from the priors. These parameters are used to create one fake dataset. Both the true underlying parameters and the simulated values are saved. Then, we create graphs showing the prior predictive distribution of the simulated values to check whether our priors fit our general expectations about the data. Next, we perform checks of computational faithfulness and model sensitivity as proposed by Schad, Betancourt and Vasishth (2020) and implemented in the SBC package. We create models for each of the simulated datasets. Last, we calculate performance metrics for each of these models, focusing on the population-level parameters.

These evaluations were performed for models with two groups; however, they should still apply to the same models with three groups.

1.4 Preparation

First, we create the data. Then, all predictors are set to sum contrasts.

# load the trial order
df.trl = read_csv('data/PAL_scheme.csv', show_col_types = F) %>%
  select(trl, emo, tone, difficulty, phase, expected)

# create the data
df.pal = data.frame()
groups = c("A", "B")
nos    = 22 # number of subjects per group
for (g in groups) {
  for (i in 1:nos) {
    df.pal = rbind(df.pal, 
                   df.trl %>% mutate(subID = sprintf("%s%02d", g, i),
                                     diagnosis = g))
  }
}

# add bogus rts
df.pal = df.pal %>% mutate(rt.cor = 1)

# drop the neutral condition for the analysis
df.pal = df.pal %>%
  filter(expected != "neutral") %>% droplevels() %>%
  mutate_if(is.character, as.factor) %>%
  ungroup()
df.var = df.pal %>%
    group_by(subID, diagnosis, phase, expected) %>% 
    summarise(
      rt.var = sd(rt.cor, na.rm = T) + 1
    ) %>%
  filter(expected != "neutral") %>% droplevels() %>%
  mutate_if(is.character, as.factor) %>%
  ungroup()

# set and print the contrasts
contrasts(df.pal$diagnosis) = contr.sum(2)
contrasts(df.pal$diagnosis)
##   [,1]
## A    1
## B   -1
contrasts(df.pal$expected) = contr.sum(2)
contrasts(df.pal$expected)
##            [,1]
## expected      1
## unexpected   -1
contrasts(df.pal$phase) = contr.sum(3)
contrasts(df.pal$phase)
##              [,1] [,2]
## postvolatile    1    0
## prevolatile     0    1
## volatile       -1   -1
contrasts(df.pal$difficulty) = contr.sum(3)
contrasts(df.pal$difficulty)
##           [,1] [,2]
## difficult    1    0
## easy         0    1
## medium      -1   -1
contrasts(df.var$diagnosis) = contr.sum(2)
contrasts(df.var$diagnosis)
##   [,1]
## A    1
## B   -1
contrasts(df.var$expected) = contr.sum(2)
contrasts(df.var$expected)
##            [,1]
## expected      1
## unexpected   -1
contrasts(df.var$phase) = contr.sum(3)
contrasts(df.var$phase)
##              [,1] [,2]
## postvolatile    1    0
## prevolatile     0    1
## volatile       -1   -1

2 Reaction time variance

In the preregistration, we noted the following population-level effects for the model of the reaction time variances: group, expectancy, phase and difficulty. However, the posterior fit for this model was suboptimal, therefore, we dropped the predictor difficulty.

2.1 Model setup

# figure out slopes for subjects
kable(head(df.var %>% count(subID, expected)))
subID expected n
A01 expected 3
A01 unexpected 3
A02 expected 3
A02 unexpected 3
A03 expected 3
A03 unexpected 3
kable(head(df.var %>% count(subID, phase)))
subID phase n
A01 postvolatile 2
A01 prevolatile 2
A01 volatile 2
A02 postvolatile 2
A02 prevolatile 2
A02 volatile 2
kable(head(df.var %>% count(subID, expected, phase)))
subID expected phase n
A01 expected postvolatile 1
A01 expected prevolatile 1
A01 expected volatile 1
A01 unexpected postvolatile 1
A01 unexpected prevolatile 1
A01 unexpected volatile 1
# code for filenames
code = "PAL-var"

# model formula
f.var = brms::bf(rt.var ~ diagnosis * expected * phase +
                   (expected + phase | subID)
                 )

# set weakly informative priors
priors = c(
  prior(normal(4.50, 0.50),  class = Intercept),
  prior(normal(0,    0.50),  class = sigma),
  prior(normal(0.15, 0.15),  class = sd),
  prior(normal(0,    0.25),  class = b)
)

2.2 Simulation-based calibration

# check if the SBC already exists
if (file.exists(file.path(cache_dir, sprintf("df_res_%s.rds", code)))) {
  # load in the results of the SBC
  df.results = readRDS(file.path(cache_dir, sprintf("df_res_%s.rds", code)))
  df.backend = readRDS(file.path(cache_dir, sprintf("df_div_%s.rds", code)))
  dat        = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
  # stimulate some data
  set.seed(2486)
  gen = SBC_generator_brms(f.var, data = df.var, prior = priors,
   thin = 50, warmup = 10000, refresh = 2000,
   generate_lp = TRUE, family = lognormal, init = 0.1)
  bck = SBC_backend_brms_from_generator(gen, chains = 4, thin = 1,
    warmup = 1500, iter = 6000)
  dat = generate_datasets(gen, nsim)
  saveRDS(dat, file.path(cache_dir, sprintf("dat_%s.rds", code)))
  # perform the SBC
  res = compute_SBC(dat, bck, keep_fits = F,
        cache_mode     = "results", 
        cache_location = file.path(cache_dir, sprintf("res_%s", code)))
  df.results = res$stats
  df.backend = res$backend_diagnostics
  saveRDS(df.results, file = file.path(cache_dir, 
                                       paste0("df_res_", code, ".rds")))
  saveRDS(df.backend, file = file.path(cache_dir, 
                                       paste0("df_div_", code, ".rds")))
}

We start by investigating the rhats and the number of divergent samples. This shows that only 0 of 250 simulations had at least one parameter that had an rhat of at least 1.05 but 14 models had divergent samples (mean = 5.93). This suggests that this model performs well, but we will check for divergence issues in the final models.

Next, we can plot the simulated values to perform prior predictive checks.

# create a matrix out of generated data
dvfakemat = matrix(NA, nrow(dat[['generated']][[1]]), length(dat[['generated']])) 
for (i in 1:length(dat[['generated']])) {
  dvfakemat[,i] = dat[['generated']][[i]][[1]]
}

# compute one histogram per simulated data-set 
dvfakematH = dvfakemat
dvfakematH[dvfakematH > 1000] = 1000
breaks = seq(0, max(dvfakematH, na.rm=T), length.out = 101) 
binwidth = breaks[2] - breaks[1]
histmat = matrix(NA, ncol = dim(dvfakemat)[2], nrow = length(breaks)-1) 
for (i in 1:dim(dvfakemat)[2]) {
  histmat[,i] = hist(dvfakematH[,i], breaks = breaks, plot = F)$counts 
}
# for each bin, compute quantiles across histograms 
probs = seq(0.1, 0.9, 0.1) 
quantmat = as.data.frame(matrix(NA, nrow=dim(histmat)[1], ncol = length(probs)))
names(quantmat) = paste0("p", probs)
for (i in 1:dim(histmat)[1]) {
  quantmat[i,] = quantile(histmat[i,], p = probs)
}
quantmat$x = breaks[2:length(breaks)] - binwidth/2 # add bin mean 
p1 = ggplot(data = quantmat, aes(x = x)) + 
  geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) + 
  geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) + 
  geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) + 
  geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) + 
  geom_line(aes(y = p0.5), colour = c_dark, linewidth = 1) + 
  labs(title = "Distribution of simulated values", y = "", x = "") +
  theme_bw()

tmpM = apply(dvfakemat, 2, mean) # mean 
tmpSD = apply(dvfakemat, 2, sd) 
p2 = ggplot() + 
  stat_bin(aes(x = tmpM), fill = c_dark)  + 
  labs(x = "Mean", title = "Means of simulated values") +
  theme_bw()
p3 = ggplot() + 
  stat_bin(aes(x = tmpSD), fill = c_dark) + 
  labs(x = "SD", title = "SDs of simulated values") +
  theme_bw()
p = ggarrange(p1, 
  ggarrange(p2, p3, ncol = 2, labels = c("B", "C")), 
  nrow = 2, labels = "A")
annotate_figure(p, top = text_grob("Prior predictive checks", face = "bold", size = 14))

Subfigure A shows the distribution of the simulated data with bluer bands being more likely than greyer bands. It shows a general distribution that fits our expectations, even though there are a few values that are unreastically large. The same applies to the distribution of the means and standard deviations in the simulated datasets in Subfigure B. We go ahead with these priors and check the results of the SBC. We only plot the results from the models that had no divergence issues.

# get simulation numbers with issues
mx_rnk = max(df.results$max_rank)
check = merge(df.results %>% 
    group_by(sim_id) %>% 
      summarise(rhat = max(rhat, na.rm = T),
                max_rank = mean(max_rank)) %>% 
    filter(rhat >= 1.05 | max_rank != mx_rnk), 
  df.backend %>% filter(n_divergent > 0), all = T)

# plot SBC with functions from the SBC package focusing on population-level parameters
a = 0.5
df.results.b = df.results %>% 
  filter(substr(variable, 1, 2) == "b_") %>% 
  filter(!(sim_id %in% check$sim_id))
p1 = plot_ecdf_diff(df.results.b) + theme_bw() + theme(legend.position = "none") +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3))
p2 = plot_rank_hist(df.results.b, bins = 20) + theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3))

p = ggarrange(p1, p2, labels = "AUTO", ncol = 1, nrow = 2)
annotate_figure(p, top = text_grob("Computational faithfulness and model sensitivity", 
                                   face = "bold", size = 14))

Next, we check the ranks of the parameters. If the model is unbiased, these should be uniformly distributed (Schad, Betancourt and Vasishth, 2020). The sample empirical cumulative distribution function (ECDF) shows deviations from the theoretical distribution (95%) and the rank histogram also shows ranks outside the 95% expected range. However, we judge this to be acceptable as we cannot identify clear bias patterns as described in the information accompanying the SBC package: https://hyunjimoon.github.io/SBC/articles/rank_visualizations.html

p3 = plot_sim_estimated(df.results.b, alpha = .8) + theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3))
p4 = plot_contraction(df.results.b, 
                      prior_sd = 
                        setNames(c(0.50, 
                                   rep(0.25,
                                       length(unique(df.results.b$variable))-1)), 
                                 unique(df.results.b$variable))) +
  theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3))

p = ggarrange(p3, p4, labels = "AUTO", ncol = 1, nrow = 2)
annotate_figure(p, top = text_grob("Computational faithfulness and model sensitivity", 
                                   face = "bold", size = 14))

Then, we investigated the relationship between the simulated true parameters and the posterior estimates. Although there are individual values diverging from the expected pattern, most parameters were recovered successfully within an uncertainty interval of alpha = 0.05.

Last, we explore the z-score and the posterior contraction of our population-level predictors. The z-score “determines the distance of the posterior mean from the true simulating parameter”, while the posterior contraction “estimates how much prior uncertainty is reduced in the posterior estimation” (Schad, Betancourt and Vasisth, 2020). The contraction reveals very narrow posterior standard deviations with slightly lower contraction scores than we would want in the optimal case, however, we judge this as acceptable.

3 Reaction times in correct trials

In the preregistration, we noted the following population-level effects for the model of the reaction time variances: group, expectancy, phase and difficulty; as well as the group-level predictores subject and trials.

3.1 Model setup

# figure out slopes for subjects
kable(head(df.pal %>% count(subID, expected)))
subID expected n
A01 expected 240
A01 unexpected 48
A02 expected 240
A02 unexpected 48
A03 expected 240
A03 unexpected 48
kable(head(df.pal %>% count(subID, phase)))
subID phase n
A01 postvolatile 72
A01 prevolatile 72
A01 volatile 144
A02 postvolatile 72
A02 prevolatile 72
A02 volatile 144
kable(head(df.pal %>% count(subID, difficulty)))
subID difficulty n
A01 difficult 96
A01 easy 96
A01 medium 96
A02 difficult 96
A02 easy 96
A02 medium 96
kable(head(df.pal %>% count(subID, expected, phase)))
subID expected phase n
A01 expected postvolatile 60
A01 expected prevolatile 60
A01 expected volatile 120
A01 unexpected postvolatile 12
A01 unexpected prevolatile 12
A01 unexpected volatile 24
kable(head(df.pal %>% count(subID, expected, difficulty)))
subID expected difficulty n
A01 expected difficult 80
A01 expected easy 80
A01 expected medium 80
A01 unexpected difficult 16
A01 unexpected easy 16
A01 unexpected medium 16
kable(head(df.pal %>% count(subID, phase, difficulty)))
subID phase difficulty n
A01 postvolatile difficult 24
A01 postvolatile easy 24
A01 postvolatile medium 24
A01 prevolatile difficult 24
A01 prevolatile easy 24
A01 prevolatile medium 24
kable(head(df.pal %>% count(subID, expected, phase, difficulty)))
subID expected phase difficulty n
A01 expected postvolatile difficult 20
A01 expected postvolatile easy 20
A01 expected postvolatile medium 20
A01 expected prevolatile difficult 20
A01 expected prevolatile easy 20
A01 expected prevolatile medium 20
code = "PAL-rt"

# set the formula
f.pal = brms::bf(rt.cor ~ diagnosis * expected * phase * difficulty +
                   (expected * phase * difficulty | subID) + (1 | trl))

# informative priors based Lawson et al. and Schad, Betancourt & Vasishth (2019)
priors = c(
  prior(normal(6.0,   0.3),   class = Intercept),
  prior(normal(0.0,   0.5),   class = sigma),
  prior(normal(0,     0.1),   class = sd),
  prior(lkj(2),               class = cor),
  prior(normal(100, 100.0),   class = ndt), 
  prior(normal(0.00,  0.04),  class = b)
)

3.2 Simulation-based calibration

Since this model is so complex and the SBC takes ages to run, we only run 125 simulations - these already took almost two weeks to run. We ran them in a separate bash script in smaller batches. Both scripts, R and bash, are in the helper folder.

if (!(file.exists(file.path(cache_dir, sprintf("df_res_%s.rds", code)))) | 
    !(file.exists(file.path(cache_dir, sprintf("df_div_%s.rds", code))))) {
  # read in the data
  dat = readRDS(file = file.path(cache_dir, paste0("dat_", code, ".rds")))
  
  # combine the SBC files to one
  ls.res = list.files(path = cache_dir, pattern = sprintf("^res_%s*", code))
  for (i in 1:length(ls.res)) {
    tmp = readRDS(file.path(cache_dir, ls.res[i]))
    if (i == 1) {
      res = tmp$result
    } else {
      res = bind_results(res, tmp$result)
    }
  }
  # put in separate variables
  df.results = res$stats
  df.backend = res$backend_diagnostics
  
  # save the output
  saveRDS(df.results, file.path(cache_dir, sprintf("df_res_%s.rds", code)))
  saveRDS(df.backend, file.path(cache_dir, sprintf("df_div_%s.rds", code)))
} else {
  # load in the results of the SBC
  df.results = readRDS(file.path(cache_dir, sprintf("df_res_%s.rds", code)))
  df.backend = readRDS(file.path(cache_dir, sprintf("df_div_%s.rds", code)))
  # read in the data
  dat = readRDS(file = file.path(cache_dir, paste0("dat_", code, ".rds")))
}

Again, we start by investigating the rhats and the number of divergent samples. This shows that 2 of 125 simulations had at least one parameter that had an rhat of at least 1.05, and 0 models had divergent samples. This suggests that this model performs well and we can continue with our checks by plotting the simulated values to perform prior predictive checks.

# create a matrix out of generated data
dvfakemat = matrix(NA, nrow(dat[['generated']][[1]]), length(dat[['generated']])) 
for (i in 1:length(dat[['generated']])) {
  dvfakemat[,i] = dat[['generated']][[i]][[1]]
}

# compute one histogram per simulated data-set 
dvfakematH = dvfakemat
dvfakematH[dvfakematH > 2000] = 2000
breaks = seq(0, max(dvfakematH, na.rm=T), length.out = 101) 
binwidth = breaks[2] - breaks[1]
histmat = matrix(NA, ncol = dim(dvfakemat)[2], nrow = length(breaks)-1) 
for (i in 1:dim(dvfakemat)[2]) {
  histmat[,i] = hist(dvfakematH[,i], breaks = breaks, plot = F)$counts 
}
# for each bin, compute quantiles across histograms 
probs = seq(0.1, 0.9, 0.1) 
quantmat = as.data.frame(matrix(NA, nrow=dim(histmat)[1], ncol = length(probs)))
names(quantmat) = paste0("p", probs)
for (i in 1:dim(histmat)[1]) {
  quantmat[i,] = quantile(histmat[i,], p = probs)
}
quantmat$x = breaks[2:length(breaks)] - binwidth/2 # add bin mean 
p1 = ggplot(data = quantmat, aes(x = x)) + 
  geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) + 
  geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) + 
  geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) + 
  geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) + 
  geom_line(aes(y = p0.5), colour = c_dark, linewidth = 1) + 
  labs(title = "Distribution of simulated values", y = "", x = "") +
  theme_bw()

tmpM = apply(dvfakemat, 2, mean) # mean 
tmpSD = apply(dvfakemat, 2, sd) 
p2 = ggplot() + 
  stat_bin(aes(x = tmpM), fill = c_dark)  + 
  labs(x = "Mean", title = "Means of simulated values") +
  theme_bw()
p3 = ggplot() + 
  stat_bin(aes(x = tmpSD), fill = c_dark) + 
  labs(x = "SD", title = "SDs of simulated values") +
  theme_bw()
p = ggarrange(p1, 
  ggarrange(p2, p3, ncol = 2, labels = c("B", "C")), 
  nrow = 2, labels = "A")
annotate_figure(p, top = text_grob("Prior predictive checks", face = "bold", size = 14))

We are happy with the distribution of the simulated values.

# get simulation numbers with issues
rank = max(df.results$max_rank)
check = merge(df.results %>% 
    group_by(sim_id) %>% 
      summarise(
        rhat = max(rhat, na.rm = T), 
        mean_rank = mean(max_rank)
        ) %>% 
    filter(rhat >= 1.05 | mean_rank != rank), 
  df.backend %>% filter(n_divergent > 0), all = T)

# plot SBC with functions from the SBC package focusing on population-level parameters
df.results.b = df.results %>% 
  filter(substr(variable, 1, 2) == "b_") %>% 
  filter(!(sim_id %in% check$sim_id))
plot_ecdf_diff(df.results.b) + theme_bw() + theme(legend.position = "none") +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: empirical cumulative distribution function")

plot_rank_hist(df.results.b, bins = 20) + theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: rank histograms")

plot_sim_estimated(df.results.b, alpha = 0.5) + theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: simulated and estimated values")

plot_contraction(df.results.b, 
                      prior_sd = setNames(
                        c(as.numeric(
                          gsub(".*, (.+)\\).*", "\\1", 
                               priors[priors$class == "Intercept",]$prior)), 
                          rep(0.04, times = (length(unique(df.results.b$variable))-1))), 
                                          unique(df.results.b$variable))) +
  theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 5)) +
  ggtitle("Computational faithfulness: contraction and z-values")

While some of the ECDFs show discrepancies from the theoretical distribution, we judge all evaluated parameters to be acceptable.

4 Pupil sizes

4.1 Model setup

# model formula
f.pup = brms::bf(rel_pupil ~ diagnosis * expected + rts + 
                   (1 | subID)
                 )

# set weakly informative priors
priors = c(
  prior(normal(0,    0.10),  class = Intercept),
  prior(normal(0.15, 0.15),  class = sigma),
  prior(normal(0.15, 0.15),  class = sd),
  prior(normal(0,    0.05),  class = b)
)

kable(priors)
prior class coef group resp dpar nlpar lb ub source
normal(0, 0.1) Intercept NA NA user
normal(0.15, 0.15) sigma NA NA user
normal(0.15, 0.15) sd NA NA user
normal(0, 0.05) b NA NA user
code = "PAL-pup-td15"

# extract one dataset from the rt.cor for this
df = dat$generated[[24]] %>%
  # add the pupil sizes
  mutate(
    rel_pupil = rnorm(nrow(dat$generated[[24]]), sd = 0.10),
    rts       = rt.cor
  ) %>% 
  # remove some of the trials to account for data loss
  filter(trl <= 252) %>%
  select(rel_pupil, diagnosis, expected, subID, rts) %>%
  mutate_if(is.character, as.factor)

# print the contrasts
contrasts(df$diagnosis) = contr.sum(2)
contrasts(df$diagnosis)
##   [,1]
## A    1
## B   -1
contrasts(df$expected)  = contr.sum(2)
contrasts(df$expected)
##            [,1]
## expected      1
## unexpected   -1

4.2 Simulation-based calibration

# check if the SBC already exists
if (file.exists(file.path(cache_dir, sprintf("df_res_%s.rds", code)))) {
  # load in the results of the SBC
  df.results = readRDS(file.path(cache_dir, sprintf("df_res_%s.rds", code)))
  df.backend = readRDS(file.path(cache_dir, sprintf("df_div_%s.rds", code)))
  dat        = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
  # stimulate some data
  set.seed(2486)
  gen = SBC_generator_brms(f.pup, data = df, prior = priors, 
   thin = 50, warmup = 10000, refresh = 2000,
   generate_lp = TRUE, init = 0.1)
  if (file.exists(file.path(cache_dir, sprintf("dat_%s.rds", code)))) {
    dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
  } else {
    dat = generate_datasets(gen, nsim)
    saveRDS(dat, file.path(cache_dir, sprintf("dat_%s.rds", code)))
  }
  # perform the SBC
  bck = SBC_backend_brms_from_generator(gen, chains = 4, thin = 1,
    warmup = 1500, iter = 6000,
    init = 0.1, 
    control = list(max_treedepth = 15))
  res = compute_SBC(dat, bck, keep_fits = F,
        cache_mode     = "results", 
        cache_location = file.path(cache_dir, sprintf("res_%s", code)))
  df.results = res$stats
  df.backend = res$backend_diagnostics
  saveRDS(df.results, file = file.path(cache_dir, 
                                       paste0("df_res_", code, ".rds")))
  saveRDS(df.backend, file = file.path(cache_dir, 
                                       paste0("df_div_", code, ".rds")))
}

Again, we start by investigating the rhats and the number of divergent samples. This shows that 134 of 250 simulations had at least one parameter that had an rhat of at least 1.05, and 0 models had divergent samples. This suggests some problems with this model. Let’s investigate this a bit further.

df.results %>% group_by(sim_id) %>% mutate(max_rhat = max(rhat, na.rm = T)) %>% filter(rhat == max_rhat & max_rhat >= 1.05) %>% group_by(variable) %>% count() %>% arrange(desc(n))
## # A tibble: 53 × 2
## # Groups:   variable [53]
##    variable                   n
##    <chr>                  <int>
##  1 Intercept                 37
##  2 b_diagnosis1              21
##  3 z_1[1,7]                   4
##  4 lprior                     3
##  5 z_1[1,14]                  3
##  6 z_1[1,17]                  3
##  7 z_1[1,6]                   3
##  8 r_subID[A04,Intercept]     2
##  9 r_subID[A07,Intercept]     2
## 10 r_subID[B01,Intercept]     2
## # ℹ 43 more rows

There are quite some large rhats for the Intercept and the diagnosis effect. We tried increasing the iterations, but it did not help. This might be connected to the maximum treedepth which is saturated in all of this simulations. We will continue for now but watch out for this in the final model.

# create a matrix out of generated data
dvfakemat = matrix(NA, nrow(dat[['generated']][[1]]), length(dat[['generated']])) 
for (i in 1:length(dat[['generated']])) {
  dvfakemat[,i] = dat[['generated']][[i]][[1]]
}

# compute one histogram per simulated data-set 
dvfakematH = dvfakemat
dvfakematH[dvfakematH > 10] = 10
dvfakematH[dvfakematH < -10] = -10
breaks = seq(min(dvfakematH, na.rm=T), max(dvfakematH, na.rm=T), length.out = 101) 
binwidth = breaks[2] - breaks[1]
histmat = matrix(NA, ncol = dim(dvfakemat)[2], nrow = length(breaks)-1) 
for (i in 1:dim(dvfakemat)[2]) {
  histmat[,i] = hist(dvfakematH[,i], breaks = breaks, plot = F)$counts 
}
# for each bin, compute quantiles across histograms 
probs = seq(0.1, 0.9, 0.1) 
quantmat = as.data.frame(matrix(NA, nrow=dim(histmat)[1], ncol = length(probs)))
names(quantmat) = paste0("p", probs)
for (i in 1:dim(histmat)[1]) {
  quantmat[i,] = quantile(histmat[i,], p = probs)
}
quantmat$x = breaks[2:length(breaks)] - binwidth/2 # add bin mean 
p1 = ggplot(data = quantmat, aes(x = x)) + 
  geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) + 
  geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) + 
  geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) + 
  geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) + 
  geom_line(aes(y = p0.5), colour = c_dark, linewidth = 1) + 
  labs(title = "Distribution of simulated values", y = "", x = "") +
  theme_bw()

tmpM = apply(dvfakemat, 2, mean) # mean 
tmpSD = apply(dvfakemat, 2, sd) 
p2 = ggplot() + 
  stat_bin(aes(x = tmpM), fill = c_dark)  + 
  labs(x = "Mean", title = "Means of simulated values") +
  theme_bw()
p3 = ggplot() + 
  stat_bin(aes(x = tmpSD), fill = c_dark) + 
  labs(x = "SD", title = "SDs of simulated values") +
  theme_bw()
p = ggarrange(p1, 
  ggarrange(p2, p3, ncol = 2, labels = c("B", "C")), 
  nrow = 2, labels = "A")
annotate_figure(p, top = text_grob("Prior predictive checks", face = "bold", size = 14))

There are quite some extreme values here that are unlikely. However, narrowing the priors also decreased the contraction. Thus, we stick with these rather wide priors.

# get simulation numbers with issues
rank = max(df.results$max_rank)
check = merge(df.results %>% 
    group_by(sim_id) %>% 
      summarise(
        rhat = max(rhat, na.rm = T), 
        mean_rank = mean(max_rank)
        ) %>% 
    filter(rhat >= 1.05 | mean_rank != rank), 
  df.backend %>% filter(n_divergent > 0), all = T)

# plot SBC with functions from the SBC package focusing on population-level parameters
df.results.b = df.results %>% 
  filter(substr(variable, 1, 2) == "b_") %>% 
  filter(!(sim_id %in% check$sim_id))
plot_ecdf_diff(df.results.b) + theme_bw() + theme(legend.position = "none") +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: empirical cumulative distribution function")

plot_rank_hist(df.results.b, bins = 20) + theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: rank histograms")

plot_sim_estimated(df.results.b, alpha = 0.5) + theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: simulated and estimated values")

plot_contraction(df.results.b, 
                      prior_sd = setNames(
                        c(as.numeric(
                          gsub(".*, (.+)\\).*", "\\1", 
                               priors[priors$class == "Intercept",]$prior)), 
                          rep(0.04, times = (length(unique(df.results.b$variable))-1))), 
                                          unique(df.results.b$variable))) +
  theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 5)) +
  ggtitle("Computational faithfulness: contraction and z-values")

Again, we have a few suboptimal values but no bias pattern. Thus, we judge this to be acceptable.

5 HGF parameters

5.1 Model setup

# model formula
f = brms::bf(p ~ diagnosis)

# set weakly informative priors
priors = c(
  prior(normal(0, 0.50),  class = Intercept),
  prior(normal(0, 0.50),  class = sigma),
  prior(normal(0, 0.25),  class = b)
)

kable(priors)
prior class coef group resp dpar nlpar lb ub source
normal(0, 0.5) Intercept NA NA user
normal(0, 0.5) sigma NA NA user
normal(0, 0.25) b NA NA user
# runs faster, so let's do some more simulations
nsim = 1000

code = "PAL-hgf"

# create dataframe
df = df.pal %>%
  select(subID, diagnosis) %>%
  distinct() %>%
  mutate(p = 1)

# print the contrasts
contrasts(df$diagnosis)
##   [,1]
## A    1
## B   -1

5.2 Simulation-based calibration

# check if the SBC already exists
if (file.exists(file.path(cache_dir, sprintf("df_res_%s.rds", code)))) {
  # load in the results of the SBC
  df.results = readRDS(file.path(cache_dir, sprintf("df_res_%s.rds", code)))
  df.backend = readRDS(file.path(cache_dir, sprintf("df_div_%s.rds", code)))
  dat        = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
  # stimulate some data
  set.seed(2486)
  gen = SBC_generator_brms(f, data = df, prior = priors, 
   thin = 50, warmup = 10000, refresh = 2000,
   generate_lp = TRUE, init = 0.1)
  if (file.exists(file.path(cache_dir, sprintf("dat_%s.rds", code)))) {
    dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
  } else {
    dat = generate_datasets(gen, nsim)
    saveRDS(dat, file.path(cache_dir, sprintf("dat_%s.rds", code)))
  }
  # perform the SBC
  bck = SBC_backend_brms_from_generator(gen, chains = 4, thin = 1,
    warmup = 1000, iter = 3000)
  res = compute_SBC(dat, bck, keep_fits = F,
        cache_mode     = "results", 
        cache_location = file.path(cache_dir, sprintf("res_%s", code)))
  df.results = res$stats
  df.backend = res$backend_diagnostics
  saveRDS(df.results, file = file.path(cache_dir, 
                                       paste0("df_res_", code, ".rds")))
  saveRDS(df.backend, file = file.path(cache_dir, 
                                       paste0("df_div_", code, ".rds")))
}

Again, we start by investigating the rhats and the number of divergent samples. This shows that 0 of 1000 simulations had at least one parameter that had an rhat of at least 1.05, and 0 models had divergent samples. This suggests no problems with this model.

# create a matrix out of generated data
dvfakemat = matrix(NA, nrow(dat[['generated']][[1]]), length(dat[['generated']])) 
for (i in 1:length(dat[['generated']])) {
  dvfakemat[,i] = dat[['generated']][[i]][[1]]
}

# compute one histogram per simulated data-set 
dvfakematH = dvfakemat
dvfakematH[dvfakematH > 10] = 10
dvfakematH[dvfakematH < -10] = -10
breaks = seq(min(dvfakematH, na.rm=T), max(dvfakematH, na.rm=T), length.out = 101) 
binwidth = breaks[2] - breaks[1]
histmat = matrix(NA, ncol = dim(dvfakemat)[2], nrow = length(breaks)-1) 
for (i in 1:dim(dvfakemat)[2]) {
  histmat[,i] = hist(dvfakematH[,i], breaks = breaks, plot = F)$counts 
}
# for each bin, compute quantiles across histograms 
probs = seq(0.1, 0.9, 0.1) 
quantmat = as.data.frame(matrix(NA, nrow=dim(histmat)[1], ncol = length(probs)))
names(quantmat) = paste0("p", probs)
for (i in 1:dim(histmat)[1]) {
  quantmat[i,] = quantile(histmat[i,], p = probs)
}
quantmat$x = breaks[2:length(breaks)] - binwidth/2 # add bin mean 
p1 = ggplot(data = quantmat, aes(x = x)) + 
  geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) + 
  geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) + 
  geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) + 
  geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) + 
  geom_line(aes(y = p0.5), colour = c_dark, linewidth = 1) + 
  labs(title = "Distribution of simulated values", y = "", x = "") +
  theme_bw()

tmpM = apply(dvfakemat, 2, mean) # mean 
tmpSD = apply(dvfakemat, 2, sd) 
p2 = ggplot() + 
  stat_bin(aes(x = tmpM), fill = c_dark)  + 
  labs(x = "Mean", title = "Means of simulated values") +
  theme_bw()
p3 = ggplot() + 
  stat_bin(aes(x = tmpSD), fill = c_dark) + 
  labs(x = "SD", title = "SDs of simulated values") +
  theme_bw()
p = ggarrange(p1, 
  ggarrange(p2, p3, ncol = 2, labels = c("B", "C")), 
  nrow = 2, labels = "A")
annotate_figure(p, top = text_grob("Prior predictive checks", face = "bold", size = 14))

We are happy with the distribution of the simulated values.

# get simulation numbers with issues
rank = max(df.results$max_rank)
check = merge(df.results %>% 
    group_by(sim_id) %>% 
      summarise(
        rhat = max(rhat, na.rm = T), 
        mean_rank = mean(max_rank)
        ) %>% 
    filter(rhat >= 1.05 | mean_rank != rank), 
  df.backend %>% filter(n_divergent > 0), all = T)

# plot SBC with functions from the SBC package focusing on population-level parameters
df.results.b = df.results %>% 
  filter(substr(variable, 1, 2) == "b_") %>% 
  filter(!(sim_id %in% check$sim_id))
plot_ecdf_diff(df.results.b) + theme_bw() + theme(legend.position = "none") +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: empirical cumulative distribution function")

plot_rank_hist(df.results.b, bins = 20) + theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: rank histograms")

plot_sim_estimated(df.results.b, alpha = 0.5) + theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: simulated and estimated values")

plot_contraction(df.results.b, 
                      prior_sd = setNames(
                        c(as.numeric(
                          gsub(".*, (.+)\\).*", "\\1", 
                               priors[priors$class == "Intercept",]$prior)), 
                          rep(0.25, times = (length(unique(df.results.b$variable))-1))), 
                                          unique(df.results.b$variable))) +
  theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 5)) +
  ggtitle("Computational faithfulness: contraction and z-values")

Again, some deviations but no bias pattern.

6 HGF learning rate update

6.1 Model setup

# model formula
f = brms::bf(p ~ diagnosis * level * change + (level + change | subID))

# set weakly informative priors
priors = c(
  prior(normal(-5,   2.00),  class = Intercept),
  prior(normal( 0.50, 0.50), class = sigma),
  prior(normal( 0,    1.00), class = b),
  prior(normal( 0.50, 0.50), class = sd),
  prior(lkj(2),              class = cor)
)

kable(priors)
prior class coef group resp dpar nlpar lb ub source
normal(-5, 2) Intercept NA NA user
normal(0.5, 0.5) sigma NA NA user
normal(0, 1) b NA NA user
normal(0.5, 0.5) sd NA NA user
lkj(2) cor NA NA user
code = "PAL-lru"

# create dataframe
df = df.pal %>%
  select(subID, diagnosis) %>%
  distinct() %>%
  mutate(p = 1)

df = rbind(df %>% mutate(level = "2", change = "pre2vol"),
           df %>% mutate(level = "2", change = "vol2post"),
           df %>% mutate(level = "3", change = "pre2vol"),
           df %>% mutate(level = "3", change = "vol2post")) %>%
  mutate_if(is.character, as.factor)

# print the contrasts
contrasts(df$diagnosis) = contr.sum(2)
contrasts(df$diagnosis)
##   [,1]
## A    1
## B   -1
contrasts(df$level) = contr.sum(2)
contrasts(df$level)
##   [,1]
## 2    1
## 3   -1
contrasts(df$change) = contr.sum(2)
contrasts(df$change)
##          [,1]
## pre2vol     1
## vol2post   -1

6.2 Simulation-based calibration

# check if the SBC already exists
if (file.exists(file.path(cache_dir, sprintf("df_res_%s.rds", code)))) {
  # load in the results of the SBC
  df.results = readRDS(file.path(cache_dir, sprintf("df_res_%s.rds", code)))
  df.backend = readRDS(file.path(cache_dir, sprintf("df_div_%s.rds", code)))
  dat        = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
  # stimulate some data
  set.seed(2486)
  gen = SBC_generator_brms(f, data = df, prior = priors, family = lognormal,
   thin = 50, warmup = 10000, refresh = 2000,
   generate_lp = TRUE, init = 0.1)
  if (file.exists(file.path(cache_dir, sprintf("dat_%s.rds", code)))) {
    dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
  } else {
    dat = generate_datasets(gen, nsim)
    saveRDS(dat, file.path(cache_dir, sprintf("dat_%s.rds", code)))
  }
  # perform the SBC
  bck = SBC_backend_brms_from_generator(gen, chains = 4, thin = 1,
    warmup = 1000, iter = 3000)
  res = compute_SBC(dat, bck, keep_fits = F,
        cache_mode     = "results", 
        cache_location = file.path(cache_dir, sprintf("res_%s", code)))
  df.results = res$stats
  df.backend = res$backend_diagnostics
  saveRDS(df.results, file = file.path(cache_dir, 
                                       paste0("df_res_", code, ".rds")))
  saveRDS(df.backend, file = file.path(cache_dir, 
                                       paste0("df_div_", code, ".rds")))
}

Again, we start by investigating the rhats and the number of divergent samples. This shows that 5 of 1000 simulations had at least one parameter that had an rhat of at least 1.05, and 21 models had divergent samples. This suggests this model performs fine.

# create a matrix out of generated data
dvfakemat = matrix(NA, nrow(dat[['generated']][[1]]), length(dat[['generated']])) 
for (i in 1:length(dat[['generated']])) {
  dvfakemat[,i] = dat[['generated']][[i]][[1]]
}

# compute one histogram per simulated data-set 
dvfakematH = dvfakemat
dvfakematH[dvfakematH > 1] = 1
#dvfakematH[dvfakematH < -10] = -10
breaks = seq(min(dvfakematH, na.rm=T), max(dvfakematH, na.rm=T), length.out = 101) 
binwidth = breaks[2] - breaks[1]
histmat = matrix(NA, ncol = dim(dvfakemat)[2], nrow = length(breaks)-1) 
for (i in 1:dim(dvfakemat)[2]) {
  histmat[,i] = hist(dvfakematH[,i], breaks = breaks, plot = F)$counts 
}
# for each bin, compute quantiles across histograms 
probs = seq(0.1, 0.9, 0.1) 
quantmat = as.data.frame(matrix(NA, nrow=dim(histmat)[1], ncol = length(probs)))
names(quantmat) = paste0("p", probs)
for (i in 1:dim(histmat)[1]) {
  quantmat[i,] = quantile(histmat[i,], p = probs)
}
quantmat$x = breaks[2:length(breaks)] - binwidth/2 # add bin mean 
p1 = ggplot(data = quantmat, aes(x = x)) + 
  geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) + 
  geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) + 
  geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) + 
  geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) + 
  geom_line(aes(y = p0.5), colour = c_dark, linewidth = 1) + 
  labs(title = "Distribution of simulated values", y = "", x = "") +
  theme_bw()

tmpM = apply(dvfakemat, 2, mean) # mean 
tmpSD = apply(dvfakemat, 2, sd) 
p2 = ggplot() + 
  stat_bin(aes(x = tmpM), fill = c_dark)  + 
  labs(x = "Mean", title = "Means of simulated values") +
  theme_bw()
p3 = ggplot() + 
  stat_bin(aes(x = tmpSD), fill = c_dark) + 
  labs(x = "SD", title = "SDs of simulated values") +
  theme_bw()
p = ggarrange(p1, 
  ggarrange(p2, p3, ncol = 2, labels = c("B", "C")), 
  nrow = 2, labels = "A")
annotate_figure(p, top = text_grob("Prior predictive checks", face = "bold", size = 14))

Again, some extreme values but we stick with the wider priors.

# get simulation numbers with issues
rank = max(df.results$max_rank)
check = merge(df.results %>% 
    group_by(sim_id) %>% 
      summarise(
        rhat = max(rhat, na.rm = T), 
        mean_rank = mean(max_rank)
        ) %>% 
    filter(rhat >= 1.05 | mean_rank != rank), 
  df.backend %>% filter(n_divergent > 0), all = T)

# plot SBC with functions from the SBC package focusing on population-level parameters
df.results.b = df.results %>% 
  filter(substr(variable, 1, 2) == "b_") %>% 
  filter(!(sim_id %in% check$sim_id))
plot_ecdf_diff(df.results.b) + theme_bw() + theme(legend.position = "none") +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: empirical cumulative distribution function")

plot_rank_hist(df.results.b, bins = 20) + theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: rank histograms")

plot_sim_estimated(df.results.b, alpha = 0.5) + theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: simulated and estimated values")

plot_contraction(df.results.b, 
                      prior_sd = setNames(
                        c(as.numeric(
                          gsub(".*, (.+)\\).*", "\\1", 
                               priors[priors$class == "Intercept",]$prior)), 
                          rep(1, times = (length(unique(df.results.b$variable))-1))), 
                                          unique(df.results.b$variable))) +
  theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 5)) +
  ggtitle("Computational faithfulness: contraction and z-values")

We judge this to be acceptable.

7 HGF Bernoulli

7.1 Model setup

# model formula
f = brms::bf(p ~ s1 + s1 + s3 + s4 + s5 + s6 + s7 )

# set weakly informative priors
priors = c(
  prior(normal(0.5, 0.50),  class = Intercept),
  prior(normal(0,   1.00), class = b)
)

kable(priors)
prior class coef group resp dpar nlpar lb ub source
normal(0.5, 0.5) Intercept NA NA user
normal(0, 1) b NA NA user
code = "PAL-ber"

# create dataframe
df = df.pal %>%
  select(subID, diagnosis) %>%
  distinct() %>%
  mutate(p  = if_else(diagnosis == "A", 1, 0),
         s1 = rnorm(length(unique(df.pal$subID))),
         s2 = rnorm(length(unique(df.pal$subID))),
         s3 = rnorm(length(unique(df.pal$subID))),
         s4 = rnorm(length(unique(df.pal$subID))),
         s5 = rnorm(length(unique(df.pal$subID))),
         s6 = rnorm(length(unique(df.pal$subID))),
         s7 = rnorm(length(unique(df.pal$subID)))
         )

7.2 Simulation-based calibration

# check if the SBC already exists
if (file.exists(file.path(cache_dir, sprintf("df_res_%s.rds", code)))) {
  # load in the results of the SBC
  df.results = readRDS(file.path(cache_dir, sprintf("df_res_%s.rds", code)))
  df.backend = readRDS(file.path(cache_dir, sprintf("df_div_%s.rds", code)))
  dat        = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
  # stimulate some data
  set.seed(2486)
  gen = SBC_generator_brms(f, data = df, prior = priors, family = bernoulli,
   thin = 50, warmup = 10000, refresh = 2000,
   generate_lp = TRUE, init = 0.1)
  if (file.exists(file.path(cache_dir, sprintf("dat_%s.rds", code)))) {
    dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
  } else {
    dat = generate_datasets(gen, nsim)
    saveRDS(dat, file.path(cache_dir, sprintf("dat_%s.rds", code)))
  }
  # perform the SBC
  bck = SBC_backend_brms_from_generator(gen, chains = 4, thin = 1,
    warmup = 1000, iter = 3000)
  res = compute_SBC(dat, bck, keep_fits = F,
        cache_mode     = "results", 
        cache_location = file.path(cache_dir, sprintf("res_%s", code)))
  df.results = res$stats
  df.backend = res$backend_diagnostics
  saveRDS(df.results, file = file.path(cache_dir, 
                                       paste0("df_res_", code, ".rds")))
  saveRDS(df.backend, file = file.path(cache_dir, 
                                       paste0("df_div_", code, ".rds")))
}

Again, we start by investigating the rhats and the number of divergent samples. This shows that 0 of 250 simulations had at least one parameter that had an rhat of at least 1.05, and 0 models had divergent samples. This suggests this model performs fine.

# create a matrix out of generated data
dvfakemat = matrix(NA, nrow(dat[['generated']][[1]]), length(dat[['generated']])) 
for (i in 1:length(dat[['generated']])) {
  dvfakemat[,i] = dat[['generated']][[i]][[1]]
}

# compute one histogram per simulated data-set 
options = c(0, 1)
histmat = matrix(NA, ncol = dim(dvfakemat)[2], length(options)) 
for (i in 1:dim(dvfakemat)[2]) {
  for (j in 1:length(options))
  {
    histmat[j,i] = sum(dvfakemat[,i] == options[j])
  }
}
# for each bin, compute quantiles across histograms 
probs = seq(0.1, 0.9, 0.1) 
quantmat= as.data.frame(matrix(NA, nrow=dim(histmat)[1], ncol = length(probs)))
names(quantmat) = paste0("p", probs)
for (i in 1:dim(histmat)[1]) {
  quantmat[i,] = quantile(histmat[i,], p = probs)
}
quantmat$x = c("group 0", "group 1")
p1 = ggplot(data = quantmat, aes(x = x)) + 
  geom_bar(aes(y = p0.9), fill = c_light, stat = "identity") + 
  geom_bar(aes(y = p0.8), fill = c_light_highlight, stat = "identity") + 
  geom_bar(aes(y = p0.7), fill = c_mid, stat = "identity") + 
  geom_bar(aes(y = p0.6), fill = c_mid_highlight, stat = "identity") + 
  geom_bar(aes(y = p0.5), fill = c_dark, stat = "identity") + 
  labs(title = "Prior predictive distribution", y = "", x = "") +
  theme_bw()

tmpM = apply(dvfakemat, 2, mean) # mean 
tmpSD = apply(dvfakemat, 2, sd) 
p2 = ggplot() + 
  stat_bin(aes(x = tmpM), fill = c_dark)  + 
  labs(x = "Mean", title = "Means of simulated values") +
  theme_bw()
p3 = ggplot() + 
  stat_bin(aes(x = tmpSD), fill = c_dark) + 
  labs(x = "SD", title = "SDs of simulated values") +
  theme_bw()
p = ggarrange(p1, 
  ggarrange(p2, p3, ncol = 2, labels = c("B", "C")), 
  nrow = 2, labels = "A")
annotate_figure(p, top = text_grob("Prior predictive checks", face = "bold", size = 14))

There is a slight bias for one group, but we judge this to be fine.

# get simulation numbers with issues
rank = max(df.results$max_rank)
check = merge(df.results %>% 
    group_by(sim_id) %>% 
      summarise(
        rhat = max(rhat, na.rm = T), 
        mean_rank = mean(max_rank)
        ) %>% 
    filter(rhat >= 1.05 | mean_rank != rank), 
  df.backend %>% filter(n_divergent > 0), all = T)

# plot SBC with functions from the SBC package focusing on population-level parameters
df.results.b = df.results %>% 
  filter(substr(variable, 1, 2) == "b_") %>% 
  filter(!(sim_id %in% check$sim_id))
plot_ecdf_diff(df.results.b) + theme_bw() + theme(legend.position = "none") +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: empirical cumulative distribution function")

plot_rank_hist(df.results.b, bins = 20) + theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: rank histograms")

plot_sim_estimated(df.results.b, alpha = 0.5) + theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: simulated and estimated values")

plot_contraction(df.results.b, 
                      prior_sd = setNames(
                        c(as.numeric(
                          gsub(".*, (.+)\\).*", "\\1", 
                               priors[priors$class == "Intercept",]$prior)), 
                          rep(1, times = (length(unique(df.results.b$variable))-1))), 
                                          unique(df.results.b$variable))) +
  theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 5)) +
  ggtitle("Computational faithfulness: contraction and z-values")

We judge this to be acceptable.

8 DDM drift rate

8.1 Model setup

# model formula
f = brms::bf(p ~ diagnosis * phase  + (1 | subID) )

# set weakly informative priors
priors = c(
  prior(normal( 2, 1.50),  class = Intercept),
  prior(normal( 0, 0.50), class = sigma),
  prior(normal( 0, 0.50), class = b),
  prior(normal( 0, 0.50), class = sd)
)

kable(priors)
prior class coef group resp dpar nlpar lb ub source
normal(2, 1.5) Intercept NA NA user
normal(0, 0.5) sigma NA NA user
normal(0, 0.5) b NA NA user
normal(0, 0.5) sd NA NA user
code = "PAL-ddm"

# create dataframe
df = df.pal %>%
  select(subID, diagnosis, phase) %>%
  distinct() %>%
  mutate(p = 1)

# print the contrasts
contrasts(df$diagnosis) = contr.sum(2)
contrasts(df$diagnosis)
##   [,1]
## A    1
## B   -1
contrasts(df$phase) = contr.sum(3)
contrasts(df$phase)
##              [,1] [,2]
## postvolatile    1    0
## prevolatile     0    1
## volatile       -1   -1

8.2 Simulation-based calibration

# check if the SBC already exists
if (file.exists(file.path(cache_dir, sprintf("df_res_%s.rds", code)))) {
  # load in the results of the SBC
  df.results = readRDS(file.path(cache_dir, sprintf("df_res_%s.rds", code)))
  df.backend = readRDS(file.path(cache_dir, sprintf("df_div_%s.rds", code)))
  dat        = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
  # stimulate some data
  set.seed(2486)
  gen = SBC_generator_brms(f, data = df, prior = priors, 
   thin = 50, warmup = 10000, refresh = 2000,
   generate_lp = TRUE, init = 0.1)
  if (file.exists(file.path(cache_dir, sprintf("dat_%s.rds", code)))) {
    dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
  } else {
    dat = generate_datasets(gen, nsim)
    saveRDS(dat, file.path(cache_dir, sprintf("dat_%s.rds", code)))
  }
  # perform the SBC
  bck = SBC_backend_brms_from_generator(gen, chains = 4, thin = 1,
    warmup = 1000, iter = 3000)
  res = compute_SBC(dat, bck, keep_fits = F,
        cache_mode     = "results", 
        cache_location = file.path(cache_dir, sprintf("res_%s", code)))
  df.results = res$stats
  df.backend = res$backend_diagnostics
  saveRDS(df.results, file = file.path(cache_dir, 
                                       paste0("df_res_", code, ".rds")))
  saveRDS(df.backend, file = file.path(cache_dir, 
                                       paste0("df_div_", code, ".rds")))
}

Again, we start by investigating the rhats and the number of divergent samples. This shows that 3 of 1000 simulations had at least one parameter that had an rhat of at least 1.05, but 103 models had divergent samples. We will keep this in mind for the final model.

# create a matrix out of generated data
dvfakemat = matrix(NA, nrow(dat[['generated']][[1]]), length(dat[['generated']])) 
for (i in 1:length(dat[['generated']])) {
  dvfakemat[,i] = dat[['generated']][[i]][[1]]
}

# compute one histogram per simulated data-set 
dvfakematH = dvfakemat
#dvfakematH[dvfakematH > 10] = 10
#dvfakematH[dvfakematH < -10] = -10
breaks = seq(min(dvfakematH, na.rm=T), max(dvfakematH, na.rm=T), length.out = 101) 
binwidth = breaks[2] - breaks[1]
histmat = matrix(NA, ncol = dim(dvfakemat)[2], nrow = length(breaks)-1) 
for (i in 1:dim(dvfakemat)[2]) {
  histmat[,i] = hist(dvfakematH[,i], breaks = breaks, plot = F)$counts 
}
# for each bin, compute quantiles across histograms 
probs = seq(0.1, 0.9, 0.1) 
quantmat = as.data.frame(matrix(NA, nrow=dim(histmat)[1], ncol = length(probs)))
names(quantmat) = paste0("p", probs)
for (i in 1:dim(histmat)[1]) {
  quantmat[i,] = quantile(histmat[i,], p = probs)
}
quantmat$x = breaks[2:length(breaks)] - binwidth/2 # add bin mean 
p1 = ggplot(data = quantmat, aes(x = x)) + 
  geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) + 
  geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) + 
  geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) + 
  geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) + 
  geom_line(aes(y = p0.5), colour = c_dark, linewidth = 1) + 
  labs(title = "Distribution of simulated values", y = "", x = "") +
  theme_bw()

tmpM = apply(dvfakemat, 2, mean) # mean 
tmpSD = apply(dvfakemat, 2, sd) 
p2 = ggplot() + 
  stat_bin(aes(x = tmpM), fill = c_dark)  + 
  labs(x = "Mean", title = "Means of simulated values") +
  theme_bw()
p3 = ggplot() + 
  stat_bin(aes(x = tmpSD), fill = c_dark) + 
  labs(x = "SD", title = "SDs of simulated values") +
  theme_bw()
p = ggarrange(p1, 
  ggarrange(p2, p3, ncol = 2, labels = c("B", "C")), 
  nrow = 2, labels = "A")
annotate_figure(p, top = text_grob("Prior predictive checks", face = "bold", size = 14))

The distributions look good.

# get simulation numbers with issues
rank = max(df.results$max_rank)
check = merge(df.results %>% 
    group_by(sim_id) %>% 
      summarise(
        rhat = max(rhat, na.rm = T), 
        mean_rank = mean(max_rank)
        ) %>% 
    filter(rhat >= 1.05 | mean_rank != rank), 
  df.backend %>% filter(n_divergent > 0), all = T)

# plot SBC with functions from the SBC package focusing on population-level parameters
df.results.b = df.results %>% 
  filter(substr(variable, 1, 2) == "b_") %>% 
  filter(!(sim_id %in% check$sim_id))
plot_ecdf_diff(df.results.b) + theme_bw() + theme(legend.position = "none") +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: empirical cumulative distribution function")

plot_rank_hist(df.results.b, bins = 20) + theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: rank histograms")

plot_sim_estimated(df.results.b, alpha = 0.5) + theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
  ggtitle("Computational faithfulness: simulated and estimated values")

plot_contraction(df.results.b, 
                      prior_sd = setNames(
                        c(as.numeric(
                          gsub(".*, (.+)\\).*", "\\1", 
                               priors[priors$class == "Intercept",]$prior)), 
                          rep(0.5, times = (length(unique(df.results.b$variable))-1))), 
                                          unique(df.results.b$variable))) +
  theme_bw() +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 5)) +
  ggtitle("Computational faithfulness: contraction and z-values")

Slight deviations, but no clear bias pattern, thus, fine.